Chapter 8 HMSC analysis
8.2 Variance partitioning
# Compute variance partitioning
varpart=computeVariancePartitioning(m)
varpart$vals %>%
as.data.frame() %>%
rownames_to_column(var="variable") %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location")))) %>%
group_by(variable) %>%
summarise(mean=mean(value)*100,sd=sd(value)*100) %>%
tt()| variable | mean | sd |
|---|---|---|
| Random: location | 37.900015 | 25.317903 |
| logseqdepth | 56.110626 | 25.796874 |
| sex | 4.937460 | 5.612719 |
| origin | 1.051899 | 1.282563 |
# Basal tree
varpart_tree <- genome_tree
#Varpart table
varpart_table <- varpart$vals %>%
as.data.frame() %>%
rownames_to_column(var="variable") %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(genome=factor(genome, levels=rev(varpart_tree$tip.label))) %>%
mutate(variable=factor(variable, levels=rev(c("origin","sex","logseqdepth","Random: location"))))
#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate(phylum=str_remove_all(phylum, "p__"))%>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% varpart_tree$tip.label) %>%
arrange(match(genome, varpart_tree$tip.label)) %>%
mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
column_to_rownames(var = "genome") %>%
dplyr::select(phylum)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate(phylum=str_remove_all(phylum, "p__"))%>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% varpart_tree$tip.label) %>%
arrange(match(genome, varpart_tree$tip.label)) %>%
dplyr::select(phylum, colors) %>%
unique() %>%
arrange(phylum) %>%
dplyr::select(colors) %>%
pull()
# Basal ggtree
varpart_tree <- varpart_tree %>%
force.ultrametric(.,method="extend") %>%
ggtree(., size = 0.3)***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultrametric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
# Add phylum colors next to the tree tips
varpart_tree <- gheatmap(varpart_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
scale_fill_manual(values=colors_alphabetic)+
labs(fill="Phylum")
#Reset fill scale to use a different colour profile in the heatmap
varpart_tree <- varpart_tree + new_scale_fill()
# Add variance stacked barplot
vertical_tree <- varpart_tree +
scale_fill_manual(values=c("#506a96","#cccccc","#be3e2b","#f6de6c"))+
geom_fruit(
data=varpart_table,
geom=geom_bar,
mapping = aes(x=value, y=genome, fill=variable, group=variable),
pwidth = 2,
offset = 0.05,
width= 1,
orientation="y",
stat="identity")+
labs(fill="Variable")
vertical_tree8.3 Posterior estimates
# Select desired support threshold
support=0.9
negsupport=1-support
# Basal tree
postestimates_tree <- genome_tree
# Posterior estimate table
post_beta <- getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(genome=factor(genome, levels=rev(postestimates_tree$tip.label))) %>%
mutate(value = case_when(
value >= support ~ "Positive",
value <= negsupport ~ "Negative",
TRUE ~ "Neutral")) %>%
mutate(value=factor(value, levels=c("Positive","Neutral","Negative"))) %>%
pivot_wider(names_from = variable, values_from = value) %>%
#select(genome,sp_vulgaris,area_semi,area_urban,sp_vulgarisxarea_semi,sp_vulgarisxarea_urban,season_spring,season_winter,sp_vulgarisxseason_spring,sp_vulgarisxseason_winter) %>%
column_to_rownames(var="genome")
#Phylums
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate(phylum=str_remove_all(phylum, "p__")) %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% postestimates_tree$tip.label) %>%
arrange(match(genome, postestimates_tree$tip.label)) %>%
mutate(phylum = factor(phylum, levels = unique(phylum))) %>%
column_to_rownames(var = "genome") %>%
dplyr::select(phylum)
colors_alphabetic <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
mutate(phylum=str_remove_all(phylum, "p__")) %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% postestimates_tree$tip.label) %>%
arrange(match(genome, postestimates_tree$tip.label)) %>%
dplyr::select(phylum, colors) %>%
unique() %>%
arrange(phylum) %>%
dplyr::select(colors) %>%
pull()
# Basal ggtree
postestimates_tree <- postestimates_tree %>%
force.ultrametric(.,method="extend") %>%
ggtree(., size = 0.3)***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultrametric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
#Add phylum colors next to the tree tips
postestimates_tree <- gheatmap(postestimates_tree, phylum_colors, offset=-0.2, width=0.1, colnames=FALSE) +
scale_fill_manual(values=colors_alphabetic)+
labs(fill="Phylum")
#Reset fill scale to use a different colour profile in the heatmap
postestimates_tree <- postestimates_tree + new_scale_fill()
# Add posterior significant heatmap
postestimates_tree <- gheatmap(postestimates_tree, post_beta, offset=0, width=0.5, colnames=TRUE, colnames_position="top",colnames_angle=90, colnames_offset_y=1, hjust=0) +
scale_fill_manual(values=c("#be3e2b","#f4f4f4","#b2b530"))+
labs(fill="Trend")
postestimates_tree +
vexpand(.25, 1) # expand top 8.3.1 Positively associated genomes
getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(trend = case_when(
value >= support ~ "Positive",
value <= negsupport ~ "Negative",
TRUE ~ "Neutral")) %>%
filter(variable=="originTame", trend=="Positive") %>%
arrange(-value) %>%
left_join(genome_metadata,by=join_by(genome==genome)) %>%
dplyr::select(genome,phylum,class,order,family,species,value) %>%
arrange(phylum, class, family)%>%
tt()| genome | phylum | class | order | family | species | value |
|---|---|---|---|---|---|---|
| bin_CaboVerde.50 | Actinomycetota | Actinomycetia | Actinomycetales | Actinomycetaceae | NA | 0.953 |
| bin_Aruba.25 | Actinomycetota | Actinomycetia | Actinomycetales | Bifidobacteriaceae | Bifidobacterium pseudocatenulatum | 0.995 |
| bin_Aruba.2 | Actinomycetota | Actinomycetia | Actinomycetales | Bifidobacteriaceae | Bifidobacterium adolescentis | 0.994 |
| bin_Aruba.6 | Actinomycetota | Actinomycetia | Actinomycetales | Bifidobacteriaceae | Bifidobacterium longum | 0.985 |
| bin_Denmark.14 | Actinomycetota | Actinomycetia | Actinomycetales | Bifidobacteriaceae | Bifidobacterium gallinarum | 0.950 |
| bin_Aruba.47 | Actinomycetota | Actinomycetia | Mycobacteriales | Mycobacteriaceae | NA | 0.964 |
| bin_Aruba.57 | Actinomycetota | Actinomycetia | Mycobacteriales | Mycobacteriaceae | Corynebacterium pyruviciproducens | 0.931 |
| bin_Aruba.51 | Actinomycetota | Coriobacteriia | Coriobacteriales | Atopobiaceae | Parolsenella sp900544655 | 0.969 |
| bin_Denmark.44 | Actinomycetota | Coriobacteriia | Coriobacteriales | Atopobiaceae | UBA7748 sp900314535 | 0.951 |
| bin_Aruba.14 | Actinomycetota | Coriobacteriia | Coriobacteriales | Coriobacteriaceae | Collinsella stercoris | 1.000 |
| bin_Denmark.4 | Actinomycetota | Coriobacteriia | Coriobacteriales | Coriobacteriaceae | Collinsella stercoris | 1.000 |
| bin_Aruba.21 | Actinomycetota | Coriobacteriia | Coriobacteriales | Coriobacteriaceae | NA | 0.999 |
| bin_Brazil.61 | Actinomycetota | Coriobacteriia | Coriobacteriales | Coriobacteriaceae | Collinsella tanakaei | 0.998 |
| bin_Aruba.39 | Actinomycetota | Coriobacteriia | Coriobacteriales | Coriobacteriaceae | Collinsella sp902362275 | 0.992 |
| bin_Spain.84 | Actinomycetota | Coriobacteriia | Coriobacteriales | Coriobacteriaceae | Collinsella sp900555555 | 0.992 |
| bin_Brazil.151 | Actinomycetota | Coriobacteriia | Coriobacteriales | Coriobacteriaceae | Collinsella intestinalis | 0.991 |
| bin_Denmark.127 | Actinomycetota | Coriobacteriia | Coriobacteriales | Coriobacteriaceae | Collinsella sp900548365 | 0.991 |
| bin_Denmark.17 | Actinomycetota | Coriobacteriia | Coriobacteriales | Coriobacteriaceae | Collinsella bouchesdurhonensis | 0.983 |
| bin_Brazil.81 | Bacillota | Bacilli | RFN20 | CAG-826 | NA | 0.907 |
| bin_Brazil.109 | Bacillota_A | Clostridia | Oscillospirales | Butyricicoccaceae | AM07-15 sp003477405 | 0.936 |
| bin_Aruba.28 | Bacillota_A | Clostridia | Oscillospirales | Oscillospiraceae | NA | 0.995 |
| bin_Malaysia.9 | Bacillota_A | Clostridia | Oscillospirales | Oscillospiraceae | Dysosmobacter welbionis | 0.986 |
| bin_Malaysia.116 | Bacillota_A | Clostridia | Oscillospirales | Oscillospiraceae | Flavonifractor plautii | 0.981 |
| bin_Malaysia.34 | Bacillota_A | Clostridia | Oscillospirales | Oscillospiraceae | Lawsonibacter sp000177015 | 0.954 |
| bin_Aruba.54 | Bacillota_A | Clostridia | Tissierellales | Peptoniphilaceae | NA | 0.959 |
| bin_CaboVerde.35 | Bacillota_A | Clostridia | Tissierellales | Peptoniphilaceae | NA | 0.930 |
| bin_Aruba.36 | Bacillota_A | Clostridia | Tissierellales | Peptoniphilaceae | Peptoniphilus_C sp902363535 | 0.913 |
| bin_Aruba.52 | Bacillota_A | Clostridia | Tissierellales | Peptoniphilaceae | NA | 0.906 |
| bin_Malaysia.26 | Bacillota_A | Clostridia | Oscillospirales | Ruminococcaceae | NA | 0.976 |
| bin_Aruba.29 | Bacillota_A | Clostridia | Oscillospirales | Ruminococcaceae | Gemmiger variabilis_C | 0.941 |
| bin_Malaysia.103 | Bacillota_C | Negativicutes | Acidaminococcales | Acidaminococcaceae | Acidaminococcus intestini | 0.942 |
| bin_Malaysia.8 | Bacillota_C | Negativicutes | Veillonellales | Dialisteraceae | Dialister sp900543165 | 0.981 |
| bin_Denmark.60 | Bacillota_C | Negativicutes | Veillonellales | Dialisteraceae | Allisonella histaminiformans | 0.938 |
| bin_Aruba.64 | Bacillota_C | Negativicutes | Veillonellales | Megasphaeraceae | Megasphaera sp000417505 | 0.970 |
| bin_Brazil.58 | Bacillota_C | Negativicutes | Veillonellales | Megasphaeraceae | Megasphaera elsdenii | 0.969 |
| bin_CaboVerde.38 | Bacillota_C | Negativicutes | Veillonellales | Megasphaeraceae | Megasphaera sp000417505 | 0.961 |
| bin_Malaysia.81 | Bacillota_C | Negativicutes | Veillonellales | Megasphaeraceae | Megasphaera stantonii | 0.943 |
| bin_Malaysia.96 | Bacillota_C | Negativicutes | Veillonellales | Megasphaeraceae | Caecibacter sp003467125 | 0.919 |
| bin_Brazil.163 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Prevotella lascolaii | 0.992 |
| bin_CaboVerde.18 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Prevotella copri | 0.981 |
| bin_Aruba.10 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Prevotella sp900544825 | 0.978 |
| bin_Brazil.5 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Prevotella sp900540415 | 0.929 |
| bin_Spain.21 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Phocaeicola coprophilus | 0.929 |
| bin_Malaysia.18 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Prevotellamassilia sp000437675 | 0.924 |
| bin_Malaysia.117 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Prevotellamassilia sp900541335 | 0.919 |
| bin_Malaysia.77 | Bacteroidota | Bacteroidia | Bacteroidales | Bacteroidaceae | Phocaeicola sp900542985 | 0.913 |
| bin_Brazil.75 | Campylobacterota | Campylobacteria | Campylobacterales | Helicobacteraceae | NA | 0.964 |
| bin_Malaysia.61 | Campylobacterota | Campylobacteria | Campylobacterales | Helicobacteraceae | Helicobacter_C labetoulli | 0.964 |
| bin_Brazil.128 | Campylobacterota | Campylobacteria | Campylobacterales | Helicobacteraceae | NA | 0.963 |
| bin_Brazil.70 | Desulfobacterota | Desulfovibrionia | Desulfovibrionales | Desulfovibrionaceae | Mailhella sp003150275 | 0.944 |
| bin_Brazil.146 | Pseudomonadota | Alphaproteobacteria | RF32 | CAG-239 | CAG-495 sp001917125 | 0.969 |
| bin_Brazil.9 | Pseudomonadota | Alphaproteobacteria | RF32 | CAG-239 | CAG-495 sp000436375 | 0.919 |
| bin_Aruba.15 | Pseudomonadota | Gammaproteobacteria | Burkholderiales | Burkholderiaceae | NA | 0.963 |
| bin_Brazil.51 | Pseudomonadota | Gammaproteobacteria | Burkholderiales | Burkholderiaceae | Sutterella wadsworthensis_A | 0.963 |
| bin_Spain.43 | Pseudomonadota | Gammaproteobacteria | Burkholderiales | Burkholderiaceae | CAG-521 sp000437635 | 0.959 |
| bin_Brazil.64 | Pseudomonadota | Gammaproteobacteria | Burkholderiales | Burkholderiaceae | NA | 0.948 |
| bin_Brazil.92 | Pseudomonadota | Gammaproteobacteria | Burkholderiales | Burkholderiaceae | CAG-521 sp000437635 | 0.948 |
| bin_Brazil.15 | Pseudomonadota | Gammaproteobacteria | Burkholderiales | Burkholderiaceae | NA | 0.930 |
| bin_CaboVerde.33 | Pseudomonadota | Gammaproteobacteria | Enterobacterales | Succinivibrionaceae | Anaerobiospirillum sp900543125 | 0.995 |
| bin_Brazil.82 | Pseudomonadota | Gammaproteobacteria | Enterobacterales | Succinivibrionaceae | Anaerobiospirillum succiniciproducens | 0.990 |
| bin_CaboVerde.55 | Pseudomonadota | Gammaproteobacteria | Enterobacterales | Succinivibrionaceae | Anaerobiospirillum_A thomasii | 0.939 |
8.3.2 Negatively associated genomes
getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
pivot_longer(!variable, names_to = "genome", values_to = "value") %>%
mutate(trend = case_when(
value >= support ~ "Positive",
value <= negsupport ~ "Negative",
TRUE ~ "Neutral")) %>%
filter(variable=="originTame", trend=="Negative") %>%
arrange(value) %>%
left_join(genome_metadata,by=join_by(genome==genome)) %>%
dplyr::select(genome,phylum,class,order,family,species,value) %>%
arrange(phylum,class,family)%>%
tt()| genome | phylum | class | order | family | species | value |
|---|---|---|---|---|---|---|
| bin_Denmark.27 | Bacillota | Bacilli | Lactobacillales | Enterococcaceae | Enterococcus_B hirae | 0.000 |
| bin_Brazil.12 | Bacillota | Bacilli | Lactobacillales | Enterococcaceae | Enterococcus faecalis | 0.001 |
| bin_Brazil.170 | Bacillota | Bacilli | Lactobacillales | Enterococcaceae | NA | 0.014 |
| bin_CaboVerde.24 | Bacillota | Bacilli | Lactobacillales | Enterococcaceae | Enterococcus_E cecorum | 0.020 |
| bin_CaboVerde.16 | Bacillota | Bacilli | Lactobacillales | Lactobacillaceae | Limosilactobacillus reuteri | 0.000 |
| bin_Denmark.56 | Bacillota | Bacilli | Lactobacillales | Lactobacillaceae | Levilactobacillus brevis | 0.000 |
| bin_Malaysia.72 | Bacillota | Bacilli | Lactobacillales | Lactobacillaceae | Limosilactobacillus reuteri | 0.000 |
| bin_CaboVerde.60 | Bacillota | Bacilli | Lactobacillales | Lactobacillaceae | Ligilactobacillus agilis | 0.001 |
| bin_Malaysia.127 | Bacillota | Bacilli | Lactobacillales | Lactobacillaceae | Ligilactobacillus agilis | 0.001 |
| bin_Denmark.61 | Bacillota | Bacilli | Lactobacillales | Lactobacillaceae | Latilactobacillus sakei | 0.003 |
| bin_Malaysia.35 | Bacillota | Bacilli | Lactobacillales | Lactobacillaceae | Ligilactobacillus animalis | 0.003 |
| bin_Malaysia.4 | Bacillota | Bacilli | Lactobacillales | Lactobacillaceae | Ligilactobacillus ruminis | 0.035 |
| bin_CaboVerde.14 | Bacillota | Bacilli | Lactobacillales | Lactobacillaceae | Lactobacillus acidophilus | 0.067 |
| bin_Aruba.18 | Bacillota | Bacilli | Lactobacillales | Streptococcaceae | Streptococcus lutetiensis | 0.001 |
| bin_Brazil.22 | Bacillota | Bacilli | Lactobacillales | Streptococcaceae | Streptococcus pasteurianus | 0.001 |
| bin_Denmark.117 | Bacillota | Bacilli | Lactobacillales | Streptococcaceae | Lactococcus garvieae | 0.001 |
| bin_Denmark.113 | Bacillota | Bacilli | Lactobacillales | Streptococcaceae | Streptococcus alactolyticus | 0.002 |
| bin_Brazil.69 | Bacillota_A | Clostridia | Oscillospirales | Acutalibacteraceae | Clostridium_A leptum | 0.097 |
| bin_Denmark.93 | Bacillota_A | Clostridia | Lachnospirales | Anaerotignaceae | Anaerotignum sp001304995 | 0.074 |
| bin_CaboVerde.3 | Bacillota_A | Clostridia | Clostridiales | Clostridiaceae | Clostridium_P perfringens | 0.000 |
| bin_Denmark.42 | Bacillota_A | Clostridia | Clostridiales | Clostridiaceae | Clostridium_P perfringens | 0.000 |
| bin_Brazil.136 | Bacillota_A | Clostridia | Clostridiales | Clostridiaceae | Clostridium thermobutyricum | 0.041 |
| bin_Brazil.3 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Faecalimonas sp900555395 | 0.026 |
| bin_Brazil.89 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Faecalimonas sp900550235 | 0.034 |
| bin_Denmark.63 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Anaerostipes hadrus | 0.040 |
| bin_Denmark.19 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Faecalimonas sp900556835 | 0.052 |
| bin_Denmark.118 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Dorea_A longicatena | 0.067 |
| bin_Malaysia.108 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | NA | 0.081 |
| bin_Brazil.113 | Bacillota_A | Clostridia | Lachnospirales | Lachnospiraceae | Anaerobutyricum sp900754855 | 0.095 |
| bin_Denmark.50 | Bacillota_A | Clostridia | Peptostreptococcales | Peptostreptococcaceae | Peptacetobacter sp900539645 | 0.053 |
| bin_Brazil.107 | Bacillota_A | Clostridia | Monoglobales_A | UBA1381 | CAG-41 sp900066215 | 0.016 |
| bin_Brazil.159 | Fusobacteriota | Fusobacteriia | Fusobacteriales | Fusobacteriaceae | Fusobacterium_B sp900541465 | 0.040 |
| bin_Brazil.140 | Fusobacteriota | Fusobacteriia | Fusobacteriales | Fusobacteriaceae | Fusobacterium_A sp900543175 | 0.046 |
| bin_Malaysia.6 | Fusobacteriota | Fusobacteriia | Fusobacteriales | Fusobacteriaceae | Fusobacterium_B sp900545035 | 0.088 |
8.3.3 Estimate - support plot
estimate <- getPostEstimate(hM=m, parName="Beta")$mean %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
filter(variable=="originTame") %>%
pivot_longer(!variable, names_to = "genome", values_to = "mean") %>%
dplyr::select(genome,mean)
support <- getPostEstimate(hM=m, parName="Beta")$support %>%
as.data.frame() %>%
mutate(variable=m$covNames) %>%
filter(variable=="originTame") %>%
pivot_longer(!variable, names_to = "genome", values_to = "support") %>%
dplyr::select(genome,support)
phylum_colors <- read_tsv("https://raw.githubusercontent.com/earthhologenome/EHI_taxonomy_colour/main/ehi_phylum_colors.tsv") %>%
right_join(genome_metadata, by=join_by(phylum == phylum)) %>%
filter(genome %in% estimate$genome) %>%
arrange(match(genome, estimate$genome)) %>%
dplyr::select(phylum, colors) %>%
unique() %>%
arrange(phylum) %>%
dplyr::select(colors) %>%
pull()Rows: 202 Columns: 2
── Column specification ─────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: "\t"
chr (2): phylum, colors
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
inner_join(estimate,support,by=join_by(genome==genome)) %>%
mutate(significance=ifelse(support >= 0.9 | support <= 0.1,1,0)) %>%
mutate(support=ifelse(mean<0,1-support,support)) %>%
left_join(genome_metadata, by = join_by(genome == genome)) %>%
mutate(phylum = ifelse(support > 0.9, phylum, NA)) %>%
ggplot(aes(x=mean,y=support,color=phylum))+
geom_point(alpha=0.7, shape=16, size=3)+
scale_color_manual(values = phylum_colors) +
geom_vline(xintercept = 0) +
xlim(c(-0.4,0.4)) +
labs(x = "Beta regression coefficient", y = "Posterior probability") +
theme_minimal()+
theme(legend.position = "none")8.4 Correlations
#Compute the residual correlation matrix
OmegaCor = computeAssociations(m)
# Refernece tree (for sorting genomes)
genome_tree_subset <- genome_tree %>%
keep.tip(., tip=m$spNames)
#Co-occurrence matrix at the animal level
supportLevel = 0.95
toPlot = ((OmegaCor[[1]]$support>supportLevel)
+ (OmegaCor[[1]]$support<(1-supportLevel))>0)*OmegaCor[[1]]$mean
matrix <- toPlot %>%
as.data.frame() %>%
rownames_to_column(var="genome1") %>%
pivot_longer(!genome1, names_to = "genome2", values_to = "cor") %>%
mutate(genome1= factor(genome1, levels=genome_tree_subset$tip.label)) %>%
mutate(genome2= factor(genome2, levels=genome_tree_subset$tip.label)) %>%
ggplot(aes(x = genome1, y = genome2, fill = cor)) +
geom_tile() +
scale_fill_gradient2(low = "#be3e2b",
mid = "#f4f4f4",
high = "#b2b530")+
theme_void()
htree <- genome_tree_subset %>%
force.ultrametric(.,method="extend") %>%
ggtree(.)***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultrametric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
***************************************************************
* Note: *
* force.ultrametric does not include a formal method to *
* ultrametricize a tree & should only be used to coerce *
* a phylogeny that fails is.ultrametric due to rounding -- *
* not as a substitute for formal rate-smoothing methods. *
***************************************************************
#create composite figure
grid.arrange(grobs = list(matrix,vtree),
layout_matrix = rbind(c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1),
c(2,1,1,1,1,1,1,1,1,1,1,1)))8.5 Predict responses
# Select modelchain of interest
load("hmsc/fit_model1_250_10.Rdata")
gradient = c("domestic","feral")
gradientlength = length(gradient)
#Treatment-specific gradient predictions
pred <- constructGradient(m,
focalVariable = "origin",
non.focalVariables = list(logseqdepth=list(1),location=list(1))) %>%
predict(m, Gradient = ., expected = TRUE) %>%
do.call(rbind,.) %>%
as.data.frame() %>%
mutate(origin=rep(gradient,1000)) %>%
pivot_longer(!origin,names_to = "genome", values_to = "value")# weights: 9 (4 variable)
initial value 101.072331
final value 91.392443
converged
8.5.0.1 Element level
elements_table <- genome_gifts %>%
to.elements(., GIFT_db) %>%
as.data.frame()
community_elements <- pred %>%
group_by(origin, genome) %>%
mutate(row_id = row_number()) %>%
pivot_wider(names_from = genome, values_from = value) %>%
ungroup() %>%
group_split(row_id) %>%
as.list() %>%
lapply(., FUN = function(x){x %>%
dplyr::select(-row_id) %>%
column_to_rownames(var = "origin") %>%
as.data.frame() %>%
exp() %>%
t() %>%
tss() %>%
to.community(elements_table,.,GIFT_db) %>%
as.data.frame() %>%
rownames_to_column(var="origin")
})
calculate_slope <- function(x) {
lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
coef(lm_fit)[2]
}
element_predictions <- map_dfc(community_elements, function(mat) {
mat %>%
column_to_rownames(var = "origin") %>%
t() %>%
as.data.frame() %>%
rowwise() %>%
mutate(slope = calculate_slope(c_across(everything()))) %>%
dplyr::select(slope) }) %>%
t() %>%
as.data.frame() %>%
set_names(colnames(community_elements[[1]])[-1]) %>%
rownames_to_column(var="iteration") %>%
pivot_longer(!iteration, names_to="trait",values_to="value") %>%
group_by(trait) %>%
summarise(mean=mean(value),
p1 = quantile(value, probs = 0.1),
p9 = quantile(value, probs = 0.9),
positive_support = sum(value > 0)/1000,
negative_support = sum(value < 0)/1000) %>%
arrange(-positive_support)8.5.1 Positive associated functions at element level
# Positively associated
unique_funct_db<- GIFT_db[c(2,4,5,6)] %>%
distinct(Code_element, .keep_all = TRUE)
element_predictions %>%
filter(mean >0) %>%
arrange(-positive_support) %>%
filter(positive_support>=0.9) %>%
left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
arrange(Domain,Function)%>%
tt()| trait | mean | p1 | p9 | positive_support | negative_support | Domain | Function | Element |
|---|---|---|---|---|---|---|---|---|
| B0103 | 0.008804968 | 6.661144e-04 | 0.016709154 | 0.920 | 0.080 | Biosynthesis | Nucleic acid biosynthesis | UDP/UTP |
| D0504 | 0.004538946 | 4.053932e-04 | 0.009749750 | 0.917 | 0.083 | Degradation | Amino acid degradation | Methionine |
| D0507 | 0.003748020 | 2.758746e-05 | 0.007466919 | 0.901 | 0.099 | Degradation | Amino acid degradation | Leucine |
| D0906 | 0.003766165 | 1.962339e-04 | 0.008167397 | 0.926 | 0.074 | Degradation | Antibiotic degradation | Fosfomycin |
| D0208 | 0.010264085 | 2.155175e-03 | 0.018260856 | 0.942 | 0.058 | Degradation | Polysaccharide degradation | Mixed-Linkage glucans |
| D0205 | 0.012642911 | 2.073360e-03 | 0.023229274 | 0.939 | 0.061 | Degradation | Polysaccharide degradation | Pectin |
element_predictions %>%
filter(mean <0) %>%
arrange(-negative_support) %>%
filter(negative_support>=0.9) %>%
left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
arrange(Domain,Function)%>%
tt()| trait | mean | p1 | p9 | positive_support | negative_support | Domain | Function | Element |
|---|---|---|---|---|---|---|---|---|
| B0219 | -0.003901767 | -0.008826791 | -2.848873e-04 | 0.043 | 0.957 | Biosynthesis | Amino acid biosynthesis | GABA |
| B0214 | -0.021208689 | -0.038785864 | -3.609382e-03 | 0.072 | 0.928 | Biosynthesis | Amino acid biosynthesis | Glutamate |
| B0204 | -0.015449214 | -0.031753343 | -6.766342e-04 | 0.092 | 0.908 | Biosynthesis | Amino acid biosynthesis | Serine |
| B0302 | -0.004719002 | -0.009783858 | -6.819933e-04 | 0.034 | 0.966 | Biosynthesis | Amino acid derivative biosynthesis | Betaine |
| B0310 | -0.012814580 | -0.023534635 | -3.361749e-03 | 0.046 | 0.954 | Biosynthesis | Amino acid derivative biosynthesis | Tryptamine |
| B0303 | -0.011870796 | -0.021137564 | -2.779849e-03 | 0.060 | 0.940 | Biosynthesis | Amino acid derivative biosynthesis | Ectoine |
| B0309 | -0.007506723 | -0.015190011 | -1.212630e-04 | 0.098 | 0.902 | Biosynthesis | Amino acid derivative biosynthesis | Putrescine |
| B0804 | -0.016705505 | -0.030010083 | -3.381213e-03 | 0.053 | 0.947 | Biosynthesis | Aromatic compound biosynthesis | Dipicolinate |
| B0603 | -0.016546534 | -0.031262013 | -2.676829e-03 | 0.069 | 0.931 | Biosynthesis | Organic anion biosynthesis | Citrate |
| B0601 | -0.009185987 | -0.017495258 | -7.435595e-04 | 0.088 | 0.912 | Biosynthesis | Organic anion biosynthesis | Succinate |
| B0401 | -0.012323164 | -0.023303496 | -7.038749e-04 | 0.085 | 0.915 | Biosynthesis | SCFA biosynthesis | Acetate |
| B0709 | -0.002087118 | -0.003665291 | -6.444408e-04 | 0.032 | 0.968 | Biosynthesis | Vitamin biosynthesis | Tocopherol/tocotorienol (E) |
| D0517 | -0.004416113 | -0.007674836 | -1.257010e-03 | 0.026 | 0.974 | Degradation | Amino acid degradation | Ornithine |
| D0508 | -0.003149520 | -0.007229581 | -1.762986e-04 | 0.079 | 0.921 | Degradation | Amino acid degradation | Lysine |
| D0903 | -0.003862111 | -0.008824495 | -2.684042e-04 | 0.043 | 0.957 | Degradation | Antibiotic degradation | Cephalosporin |
| D0908 | -0.015343133 | -0.027690198 | -4.284755e-03 | 0.052 | 0.948 | Degradation | Antibiotic degradation | Macrolide |
| D0601 | -0.009191565 | -0.017397410 | -2.465709e-03 | 0.037 | 0.963 | Degradation | Nitrogen compound degradation | Nitrate |
| D0610 | -0.003027935 | -0.004879338 | -9.417130e-04 | 0.039 | 0.961 | Degradation | Nitrogen compound degradation | Methylamine |
| D0611 | -0.003862111 | -0.008824495 | -2.684042e-04 | 0.043 | 0.957 | Degradation | Nitrogen compound degradation | Phenylethylamine |
| D0603 | -0.001850442 | -0.003639849 | -3.827215e-04 | 0.048 | 0.952 | Degradation | Nitrogen compound degradation | Urate |
| D0606 | -0.005646982 | -0.010629987 | -1.125038e-03 | 0.060 | 0.940 | Degradation | Nitrogen compound degradation | Allantoin |
| D0612 | -0.001583928 | -0.002951867 | -1.645286e-05 | 0.098 | 0.902 | Degradation | Nitrogen compound degradation | Hypotaurine |
| D0801 | -0.001879207 | -0.002299701 | -1.080150e-04 | 0.004 | 0.996 | Degradation | Xenobiotic degradation | Toluene |
| D0802 | -0.001879207 | -0.002299701 | -1.080150e-04 | 0.004 | 0.996 | Degradation | Xenobiotic degradation | Xylene |
| D0807 | -0.004059376 | -0.008574328 | -6.102883e-04 | 0.055 | 0.945 | Degradation | Xenobiotic degradation | Catechol |
| D0817 | -0.004725492 | -0.010236141 | -6.085526e-04 | 0.056 | 0.944 | Degradation | Xenobiotic degradation | Trans-cinnamate |
| D0816 | -0.005551332 | -0.011507872 | -4.487159e-04 | 0.083 | 0.917 | Degradation | Xenobiotic degradation | Phenylacetate |
8.5.2 Negatively associated functions at element level
element_predictions %>%
filter(mean <0) %>%
arrange(-negative_support) %>%
filter(negative_support>=0.9) %>%
left_join(unique_funct_db, by = join_by(trait == Code_element))%>%
arrange(Domain,Function)%>%
tt()| trait | mean | p1 | p9 | positive_support | negative_support | Domain | Function | Element |
|---|---|---|---|---|---|---|---|---|
| B0219 | -0.003901767 | -0.008826791 | -2.848873e-04 | 0.043 | 0.957 | Biosynthesis | Amino acid biosynthesis | GABA |
| B0214 | -0.021208689 | -0.038785864 | -3.609382e-03 | 0.072 | 0.928 | Biosynthesis | Amino acid biosynthesis | Glutamate |
| B0204 | -0.015449214 | -0.031753343 | -6.766342e-04 | 0.092 | 0.908 | Biosynthesis | Amino acid biosynthesis | Serine |
| B0302 | -0.004719002 | -0.009783858 | -6.819933e-04 | 0.034 | 0.966 | Biosynthesis | Amino acid derivative biosynthesis | Betaine |
| B0310 | -0.012814580 | -0.023534635 | -3.361749e-03 | 0.046 | 0.954 | Biosynthesis | Amino acid derivative biosynthesis | Tryptamine |
| B0303 | -0.011870796 | -0.021137564 | -2.779849e-03 | 0.060 | 0.940 | Biosynthesis | Amino acid derivative biosynthesis | Ectoine |
| B0309 | -0.007506723 | -0.015190011 | -1.212630e-04 | 0.098 | 0.902 | Biosynthesis | Amino acid derivative biosynthesis | Putrescine |
| B0804 | -0.016705505 | -0.030010083 | -3.381213e-03 | 0.053 | 0.947 | Biosynthesis | Aromatic compound biosynthesis | Dipicolinate |
| B0603 | -0.016546534 | -0.031262013 | -2.676829e-03 | 0.069 | 0.931 | Biosynthesis | Organic anion biosynthesis | Citrate |
| B0601 | -0.009185987 | -0.017495258 | -7.435595e-04 | 0.088 | 0.912 | Biosynthesis | Organic anion biosynthesis | Succinate |
| B0401 | -0.012323164 | -0.023303496 | -7.038749e-04 | 0.085 | 0.915 | Biosynthesis | SCFA biosynthesis | Acetate |
| B0709 | -0.002087118 | -0.003665291 | -6.444408e-04 | 0.032 | 0.968 | Biosynthesis | Vitamin biosynthesis | Tocopherol/tocotorienol (E) |
| D0517 | -0.004416113 | -0.007674836 | -1.257010e-03 | 0.026 | 0.974 | Degradation | Amino acid degradation | Ornithine |
| D0508 | -0.003149520 | -0.007229581 | -1.762986e-04 | 0.079 | 0.921 | Degradation | Amino acid degradation | Lysine |
| D0903 | -0.003862111 | -0.008824495 | -2.684042e-04 | 0.043 | 0.957 | Degradation | Antibiotic degradation | Cephalosporin |
| D0908 | -0.015343133 | -0.027690198 | -4.284755e-03 | 0.052 | 0.948 | Degradation | Antibiotic degradation | Macrolide |
| D0601 | -0.009191565 | -0.017397410 | -2.465709e-03 | 0.037 | 0.963 | Degradation | Nitrogen compound degradation | Nitrate |
| D0610 | -0.003027935 | -0.004879338 | -9.417130e-04 | 0.039 | 0.961 | Degradation | Nitrogen compound degradation | Methylamine |
| D0611 | -0.003862111 | -0.008824495 | -2.684042e-04 | 0.043 | 0.957 | Degradation | Nitrogen compound degradation | Phenylethylamine |
| D0603 | -0.001850442 | -0.003639849 | -3.827215e-04 | 0.048 | 0.952 | Degradation | Nitrogen compound degradation | Urate |
| D0606 | -0.005646982 | -0.010629987 | -1.125038e-03 | 0.060 | 0.940 | Degradation | Nitrogen compound degradation | Allantoin |
| D0612 | -0.001583928 | -0.002951867 | -1.645286e-05 | 0.098 | 0.902 | Degradation | Nitrogen compound degradation | Hypotaurine |
| D0801 | -0.001879207 | -0.002299701 | -1.080150e-04 | 0.004 | 0.996 | Degradation | Xenobiotic degradation | Toluene |
| D0802 | -0.001879207 | -0.002299701 | -1.080150e-04 | 0.004 | 0.996 | Degradation | Xenobiotic degradation | Xylene |
| D0807 | -0.004059376 | -0.008574328 | -6.102883e-04 | 0.055 | 0.945 | Degradation | Xenobiotic degradation | Catechol |
| D0817 | -0.004725492 | -0.010236141 | -6.085526e-04 | 0.056 | 0.944 | Degradation | Xenobiotic degradation | Trans-cinnamate |
| D0816 | -0.005551332 | -0.011507872 | -4.487159e-04 | 0.083 | 0.917 | Degradation | Xenobiotic degradation | Phenylacetate |
positive <- element_predictions %>%
filter(mean >0) %>%
arrange(mean) %>%
filter(positive_support>=0.9) %>%
dplyr::select(-negative_support) %>%
rename(support=positive_support)
negative <- element_predictions %>%
filter(mean <0) %>%
arrange(mean) %>%
filter(negative_support>=0.9) %>%
dplyr::select(-positive_support) %>%
rename(support=negative_support)
bind_rows(positive,negative) %>%
left_join(GIFT_db,by=join_by(trait==Code_element)) %>%
mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
geom_point() +
geom_errorbar() +
xlim(c(-0.04,0.04)) +
geom_vline(xintercept=0) +
scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
theme_minimal() +
labs(x="Regression coefficient",y="Functional trait")table <- bind_rows(positive,negative) %>%
left_join(unique_funct_db,by=join_by(trait==Code_element)) %>%
mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait))))
table %>%
mutate(Element=factor(Element,levels=c(table$Element))) %>%
ggplot(aes(x=mean, y=fct_rev(Element), xmin=p1, xmax=p9, color=Function)) +
geom_point() +
geom_errorbar() +
xlim(c(-0.04,0.04)) +
geom_vline(xintercept=0) +
scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
theme_minimal() +
labs(x="Regression coefficient",y="Functional trait")8.5.2.1 Function level
functions_table <- elements_table %>%
to.functions(., GIFT_db) %>%
as.data.frame()
community_functions <- pred %>%
group_by(origin, genome) %>%
mutate(row_id = row_number()) %>%
pivot_wider(names_from = genome, values_from = value) %>%
ungroup() %>%
group_split(row_id) %>%
as.list() %>%
lapply(., FUN = function(x){x %>%
dplyr::select(-row_id) %>%
column_to_rownames(var = "origin") %>%
as.data.frame() %>%
exp() %>%
t() %>%
tss() %>%
to.community(functions_table,.,GIFT_db) %>%
as.data.frame() %>%
rownames_to_column(var="origin")
})#max-min option
calculate_slope <- function(x) {
lm_fit <- lm(unlist(x) ~ seq_along(unlist(x)))
coef(lm_fit)[2]
}
function_predictions <- map_dfc(community_functions, function(mat) {
mat %>%
column_to_rownames(var = "origin") %>%
t() %>%
as.data.frame() %>%
rowwise() %>%
mutate(slope = calculate_slope(c_across(everything()))) %>%
dplyr::select(slope) }) %>%
t() %>%
as.data.frame() %>%
set_names(colnames(community_functions[[1]])[-1]) %>%
rownames_to_column(var="iteration") %>%
pivot_longer(!iteration, names_to="trait",values_to="value") %>%
group_by(trait) %>%
summarise(mean=mean(value),
p1 = quantile(value, probs = 0.1),
p9 = quantile(value, probs = 0.9),
positive_support = sum(value > 0)/1000,
negative_support = sum(value < 0)/1000) %>%
arrange(-positive_support)
# Positively associated
function_predictions %>%
filter(mean >0) %>%
arrange(-positive_support) %>%
tt()| trait | mean | p1 | p9 | positive_support | negative_support |
|---|---|---|---|---|---|
| D02 | 8.380877e-03 | -0.0032042156 | 0.021012695 | 0.829 | 0.171 |
| B08 | 8.092933e-03 | -0.0036562752 | 0.019064905 | 0.803 | 0.197 |
| B01 | 1.198216e-03 | -0.0057372614 | 0.008344889 | 0.652 | 0.348 |
| S01 | 9.353796e-04 | -0.0114173781 | 0.013508236 | 0.586 | 0.414 |
| B09 | 4.002654e-05 | -0.0005537729 | 0.000499940 | 0.361 | 0.639 |
# Negatively associated
function_predictions %>%
filter(mean <0) %>%
arrange(-negative_support) %>%
tt()| trait | mean | p1 | p9 | positive_support | negative_support |
|---|---|---|---|---|---|
| D08 | -0.0011456784 | -0.002162438 | -0.0002132115 | 0.039 | 0.961 |
| B03 | -0.0106204008 | -0.018053351 | -0.0033897581 | 0.050 | 0.950 |
| D06 | -0.0029947487 | -0.006665266 | 0.0001147733 | 0.107 | 0.893 |
| B04 | -0.0084917327 | -0.018498150 | 0.0007309851 | 0.130 | 0.870 |
| B06 | -0.0067968210 | -0.016392244 | 0.0029246318 | 0.177 | 0.823 |
| D07 | -0.0121698848 | -0.031193068 | 0.0045376311 | 0.192 | 0.808 |
| D03 | -0.0041778567 | -0.013047116 | 0.0032739106 | 0.207 | 0.793 |
| D05 | -0.0014533798 | -0.007604949 | 0.0034585469 | 0.224 | 0.776 |
| S03 | -0.0094794343 | -0.032363817 | 0.0155818115 | 0.240 | 0.760 |
| B02 | -0.0034965044 | -0.012889226 | 0.0053410419 | 0.278 | 0.722 |
| D09 | -0.0019870507 | -0.008375102 | 0.0042295505 | 0.302 | 0.698 |
| S02 | -0.0046167761 | -0.014546536 | 0.0029897655 | 0.332 | 0.668 |
| B07 | -0.0034403466 | -0.016055313 | 0.0096265812 | 0.348 | 0.652 |
| D01 | -0.0001590696 | -0.005078656 | 0.0040657195 | 0.434 | 0.566 |
| B10 | -0.0000088381 | -0.000324572 | 0.0002401476 | 0.494 | 0.506 |
positive <- function_predictions %>%
filter(mean >0) %>%
arrange(mean) %>%
filter(positive_support>=0.9) %>%
dplyr::select(-negative_support) %>%
rename(support=positive_support)
negative <- function_predictions %>%
filter(mean <0) %>%
arrange(mean) %>%
filter(negative_support>=0.9) %>%
dplyr::select(-positive_support) %>%
rename(support=negative_support)
bind_rows(positive,negative) %>%
left_join(GIFT_db,by=join_by(trait==Code_function)) %>%
mutate(trait=factor(trait,levels=c(rev(positive$trait),rev(negative$trait)))) %>%
ggplot(aes(x=mean, y=fct_rev(trait), xmin=p1, xmax=p9, color=Function)) +
geom_point() +
geom_errorbar() +
xlim(c(-0.02,0.02)) +
geom_vline(xintercept=0) +
scale_color_manual(values = c("#debc14","#440526","#dc7c17","#172742","#debc14","#440526","#dc7c17","#172742","#357379","#6c7e2c","#d8dc69","#774d35","#db717d")) +
theme_minimal() +
labs(x="Regression coefficient",y="Functional trait")